Author(s): Olli Väinölä olli.vainola@student.oulu.fi
In this notebook is an small tutorial, how to create a Von Mises material without any hardening. Equations are formulated into rate depended form.
In [1]:
# imports
using ForwardDiff
using NLsolve
using PyPlot
Let's create a isotropic Hooke material.
In [2]:
"""
Create a isotropic Hooke material matrix C
More information: # http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm
Parameters
----------
E: Float
Elastic modulus
ν: Float
Poisson constant
Returns
-------
Array{Float64, (6,6)}
"""
function hookeStiffnessTensor(E, ν)
a = 1 - ν
b = 1 - 2*ν
c = 1 + ν
multiplier = E / (b * c)
return Float64[a ν ν 0 0 0;
ν a ν 0 0 0;
ν ν a 0 0 0;
0 0 0 b 0 0;
0 0 0 0 b 0;
0 0 0 0 0 b].*multiplier
end
# Pick material values
E = 200.0e3
ν = 0.3
C = hookeStiffnessTensor(E, ν)
Out[2]:
In [3]:
# using vectors with double contradiction
# http://www-2.unipv.it/compmech/teaching/available/const_mod/const_mod_mat-review_notation.pdf
M = [1 0 0 0 0 0;
0 1 0 0 0 0;
0 0 1 0 0 0;
0 0 0 2 0 0;
0 0 0 0 2 0;
0 0 0 0 0 2;]
"""
Equivalent tensile stress.
More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion
Section: Reduced von Mises equation for different stress conditions
Parameters
----------
σ: Array{Float64, 6}
Stress in Voigt notation
Returns
-------
Float
"""
function σₑ(σ)
s = σ[1:6] - 1/3 * sum([σ[1], σ[2], σ[3]]) * [1 1 1 0 0 0]'
return sqrt(3/2 * s' * M * s)[1]
end
"""
Von Mises Yield criterion
More info can be found from: http://csm.mech.utah.edu/content/wp-content/uploads/2011/10/9tutorialOnJ2Plasticity.pdf
Parameters
----------
σ: Array{Float64, 6}
Stress in Voigt notation
k: Float64
Material constant, Yield limit
Returns
-------
Float
"""
function vonMisesYield(σ, k)
σₑ(σ) - k
end
"""
Function for NLsolve. Inside this function are the equations which we want to find root.
Ψ is the yield function below. Functions defined here:
dσ - C (dϵ - dλ*dΨ/dσ) = 0
σₑ(σ) - k = 0
Parameters
----------
params: Array{Float64, 7}
Array containing values from solver
dϵ: Array{Float64, 6}
Strain rate vector in Voigt notation
C: Array{Float64, (6, 6)}
Material tensor
k: Float
Material constant, yield limit
Δt: Float
time increment
σ_begin:Array{Float64, 6}
Stress vector in Voigt notation
Returns
-------
Array{Float64, 7}, return values for solver
"""
function G(params, dϵ, C, k, σ_begin)
# Creating wrapper for gradient
yield(pars) = vonMisesYield(pars, k)
dfdσ = ForwardDiff.gradient(yield)
# Stress rate
dσ = params[1:6]
σ_tot = [vec(σ_begin); 0.0] + params
# Calculating plastic strain rate
dϵp = params[end] * dfdσ(σ_tot)
# Calculating equations
function_1 = dσ - C * (dϵ - dϵp[1:6])
function_2 = yield(σ_tot)
[vec(function_1); function_2]
end
"""
Function which calculates the stress. Also handles if any yielding happens
Parameters
----------
dϵ: Array{Float64, 6}
Strain rate vector in Voigt notation
Δt: Float
time increment
σ: Array{Float64, 6}
Last stress vector in Voigt notation
C: Array{Float64, (6, 6)}
Material tensor
k: Float
Material constant, yield limit
Returns
-------
Tuple
Plastic strain rate dϵᵖ and new stress vector σ
"""
function calculate_stress(dϵ, σ, C, k)
# Test stress
σ_tria = σ + C * dϵ
# Calculating yield
yield = vonMisesYield(σ_tria, k)
if yield > 0
# Yielding happened
# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
initial_guess = [vec(σ_tria - σ); 0.1]
f(σ_) = G(σ_, dϵ, C, k, σ)
df = ForwardDiff.jacobian(f)
# Calculating root
result = nlsolve(not_in_place(f, df), initial_guess).zero
σ[:] += result[1:6]
else
σ = σ_tria
end
return σ
end
Out[3]:
In [4]:
steps = 300
strain_max = 0.003
num_cycles = 3
ϵ_tot = zeros(Float64, (steps, 6))
ϵ_tot2 = zeros(Float64, (steps, 6))
ϵ_tot3 = zeros(Float64, (steps, 6))
# Adding only strain in x-axis and counting for the poisson effect
ϵ_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))
ϵ_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
ϵ_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
Out[4]:
In [5]:
ϵ_last = zeros(Float64, (6))
ϵᵖ = zeros(Float64, (6))
σ = zeros(Float64, (6, 1))
σy = 200.0
ss = Float64[]
ee = Float64[]
for i=1:steps
dϵ = reshape(ϵ_tot[i, :, :], (6, 1)) - ϵ_last
σ = calculate_stress(dϵ, σ, C, σy)
ϵ_last += dϵ
push!(ss, σ[1])
push!(ee, ϵ_last[1])
end
PyPlot.plot(ee, ss)
PyPlot.title("Stress-Strain curve")
PyPlot.xlabel("Strain")
PyPlot.ylabel("Stress")
Out[5]:
In [ ]: